function sb = solve_backward(param,fp,fp_parfor)

alpha = param(1:9);
TFP   = param(10:11);     
delta = param(12:14);
gamma = param(15:19);
%Utility cost of being authoritarian by neighorhood income;
delta(2) = delta(2) + param(20)*(fp.mean_income_by_neighborhoods(fp_parfor.j) - fp.mean_income_by_neighborhoods(4) );

%Discount Factor;
beta = 0.95;

%Pre-define cells used to store results;
estim_policy0=cell(fp.periods-1,1);
estim_policy1=cell(fp.periods-1,1);
            
estim_value  =cell(fp.periods-1,1);
estim_value0 =cell(fp.periods-1,1);
estim_value1 =cell(fp.periods-1,1);

%Grid points;
grid_peers=fp_parfor.grid_peers;
inv0 = fp_parfor.inv0;
inv1 = fp_parfor.inv1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Create Grid;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

hk = fp_parfor.grid ;

%Guess Distribution of Peers;

if fp_parfor.i==1
%Initial guess in the solution algorithm;    
Hc_tp1_guess{1} = max(fp_parfor.distrib_skills(:,2)).*ones( length(inv1{1}), length(hk{1}) , length(grid_peers{1})  );
Hc_tp1_guess{2} = max(fp_parfor.distrib_skills(:,3)).*ones( length(inv1{2}), length(hk{2}) , length(grid_peers{2})  );
    
else
% Previous distribution in the fixed-point algorithm;    
Hc_tp1_guess = fp_parfor.Hc_tp1_guess_prime ;
end



for t=fp.periods-1:-1:1

    grid_p=grid_peers{t};
    grid_p = sort(grid_p);


    if t==fp.periods-1     %Solve last period


        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
        %Solving the parent's problem;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

        grid_inv0 = inv0{t};
        grid_inv1 = inv1{t};

        [I0, hc0, H0] =ndgrid(grid_inv0, hk{t}, grid_p );
        [I1, hc1, H1] =ndgrid(grid_inv1, hk{t}, grid_p );

        log_Skills_tp1_0 = log(technology(t,alpha, hc0, I0, 0, H0, TFP) ) ;
        log_Skills_tp1_1 = log(technology(t,alpha, hc1, I1, 1, H1, TFP) ) ;

        %Parent's problem in the last period (equation 3) ;
        U0 = delta(1).*log(hc0).*(1 - delta(3).*0) + 1.*log(1-I0)  -delta(2).*0  + beta.*delta(1).*log_Skills_tp1_0;
        U1 = delta(1).*log(hc1).*(1 - delta(3).*1) + 1.*log(1-I1)  -delta(2).*1  + beta.*delta(1).*log_Skills_tp1_1;

        %Counterfactual when we halve the cost of parental investments;
        if fp.do_counter_type == 8
            U0 = delta(1).*log(hc0).*(1 - delta(3).*0) + 0.5.*log(1-I0)  -delta(2).*0  + beta.*delta(1).*log_Skills_tp1_0;
            U1 = delta(1).*log(hc1).*(1 - delta(3).*1) + 0.5.*log(1-I1)  -delta(2).*1  + beta.*delta(1).*log_Skills_tp1_1;
        end

        %Find the optimal parental investments by parenting style;
        [Value0,ind0] = max(U0,[],1);
        [Value1,ind1] = max(U1,[],1);


        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
        %Computing Value e Policy Functions;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
        %Value Function;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

        Value0_actual(:,:,:,t) = Value0;
        Value1_actual(:,:,:,t) = Value1;        

        %Collect State Variables;
        hc_s = hc0(1,:,:);
        H_s  = H1(1,:,:);
        State_var = [ hc_s(:), H_s(:) ];

        %Value for P=0;
        Y0 = Value0(:);
        estim_value0{t} = compute_value_function(Y0, State_var);
        %Value for P=1;
        Y1 = Value1(:);
        estim_value1{t} = compute_value_function(Y1, State_var);
        %Expected Value over preference shocks;
        Y = log( exp(Y0) + exp(Y1));
        estim_value{t}  = compute_value_function(Y, State_var);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %Policy Function (Investments)
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        %Store the optimal choice for each state (by P);
        Inv0 = I0(ind0);
        Inv1 = I1(ind1);

        Y0 = Inv0(:);
        estim_policy0{t} = compute_policy_function(Y0, State_var);
        Y1 = Inv1(:);
        estim_policy1{t} = compute_policy_function(Y1, State_var);


    else %Before Last period;

        grid_inv0 = inv0{t};
        grid_inv1 = inv1{t};

        [I0, hc0, H0] =ndgrid(grid_inv0, hk{t}, grid_p );
        [I1, hc1, H1] =ndgrid(grid_inv1, hk{t}, grid_p );

        %Distribution coming from the Solve_Forward;
        set_kids2 = fp_parfor.distrib_skills(:,t+1);
        set_PS2   = fp_parfor.distrib_PS(:,t);

        %Iteration
        set_kids = Hc_tp1_guess{t};
        set_kids1 = set_kids(:);

        %Random graph generation;
        [H1_tp1, H2_tp1] = ndgrid(set_kids1, set_kids2);

        set_PS_0 = zeros(length(set_kids1),1);
        set_PS_1 = ones(length(set_kids1),1);

        [PS1_0, PS2_0] = ndgrid(set_PS_0, set_PS2);
        [PS1_1, PS2_1] = ndgrid(set_PS_1, set_PS2);

        m_p0 = mean_peers(H1_tp1,H2_tp1,PS1_0,PS2_0,gamma,fp,set_kids2,fp_parfor,t);
        H_tp1_0 = m_p0.H_prime ;

        m_p1 = mean_peers(H1_tp1,H2_tp1,PS1_1,PS2_1,gamma,fp,set_kids2,fp_parfor,t);
        H_tp1_1 = m_p1.H_prime ;

        %Evolution of peer qualily by parenting style;
        H_prime_0 = reshape(nanmean(log(H_tp1_0),2),[size(hc0,1),size(hc0,2),size(hc0,3)]);
        H_prime_0(isnan(H_prime_0)) = nanmean(H_prime_0(:)); %Not relevant, just safety: in case peers are not defined for a state;
        H_prime_1 = reshape(nanmean(log(H_tp1_1),2),[size(hc1,1),size(hc1,2),size(hc1,3)]);
        H_prime_1(isnan(H_prime_1)) = nanmean(H_prime_1(:)); %Not relevant, just safety: in case peers are not defined for a state;

        %Evolution of a child's skills by parenting style;
        log_Skills_tp1_0 = log(technology(t,alpha, hc0, I0, 0, H0, TFP) ) ;
        log_Skills_tp1_1 = log(technology(t,alpha, hc1, I1, 1, H1, TFP) ) ;

        %Calculate Altruism part of the parents' utility (based on equation 7);

        %Calculate Altruism based on Expected Value of Child's Utility (P=0);
        temp_utility_child0 = nanmean( m_p0.U_links1,3);
        utility_child0 = nansum( temp_utility_child0,2);
        Altruism_0 = reshape( utility_child0 ,[size(hc0,1),size(hc0,2),size(hc0,3)]);

        %Calculate Altruism based on Expected Value of Child's Utility (P=1);        
        temp_utility_child1 = nanmean( m_p1.U_links1,3);
        utility_child1 = nansum( temp_utility_child1,2);
        Altruism_1 = reshape( utility_child1 ,[size(hc1,1),size(hc1,2),size(hc1,3)]);

        %Simulated Distribution of Peer Skills;
        distrib_mc_H_peers0{t} = H_tp1_0 ;
        distrib_mc_H_peers1{t} = H_tp1_1 ;


        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
        %Solving the parent's problem;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

        %Continuation value (by parenting style);
        Cont_value0 = Interp_function(estim_value{t+1},log_Skills_tp1_0,H_prime_0);
        Cont_value1 = Interp_function(estim_value{t+1},log_Skills_tp1_1,H_prime_1); 

        %Bellman Equation by Parenting Style (Equation 2 of the Paper);
        U0 = fp.lambda.*delta(1).*log(hc0).*(1 - delta(3).*0) + (1-fp.lambda).*Altruism_0 + 1.*log(1-I0)  -delta(2).*0  + beta.*Cont_value0;
        U1 = fp.lambda.*delta(1).*log(hc1).*(1 - delta(3).*1) + (1-fp.lambda).*Altruism_1 + 1.*log(1-I1)  -delta(2).*1  + beta.*Cont_value1;

        %Counterfactual when we halve the cost of parental investments;
        if fp.do_counter_type == 8
            U0 = fp.lambda.*delta(1).*log(hc0).*(1 - delta(3).*0) + (1-fp.lambda).*Altruism_0 + 0.5.*log(1-I0)  -delta(2).*0  + beta.*Cont_value0;
            U1 = fp.lambda.*delta(1).*log(hc1).*(1 - delta(3).*1) + (1-fp.lambda).*Altruism_1 + 0.5.*log(1-I1)  -delta(2).*1  + beta.*Cont_value1;
        end
        %Find the optimal parental investments by parenting style;
        [Value0,ind0] = max(U0,[],1);
        [Value1,ind1] = max(U1,[],1);

        temp_value0 = repmat( Value0 , [size(I0,1), 1, 1 ] );
        temp_value1 = repmat( Value1 , [size(I1,1), 1, 1 ] );

        %Optimal behavior about parenting style;
        I_guess0 = repmat(I0(ind0), [size(I0,1), 1, 1 ] );
        I_guess1 = repmat(I1(ind1), [size(I1,1), 1, 1 ] );

        %Logic Probability of parenting style;
        Prob_PS = (1./(1+1./exp(temp_value1 -temp_value0)));

        %Updating the guess for next-period peers skill aggregate distribution according to optimal behavior of parents;
        Hc_tp1_guess{t} = (1-Prob_PS).*technology(t,alpha, hc0, I_guess0, 0, H0, TFP) + Prob_PS.*technology(t,alpha, hc1, I_guess1, 1, H1, TFP);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
        %Computing Value e Policy Functions;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
        %Value Function;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

        Value0_actual(:,:,:,t) = Value0;
        Value1_actual(:,:,:,t) = Value1;        

        %Collect State Variables;
        hc_s = hc1(1,:,:);
        H_s  = H1(1,:,:);
        State_var = [ hc_s(:), H_s(:) ];

        %Value for P=0;
        Y0 = Value0(:);
        estim_value0{t} = compute_value_function(Y0, State_var);
        %Value for P=1;
        Y1 = Value1(:);
        estim_value1{t} = compute_value_function(Y1, State_var);

        %Expected Value over preference shocks;
        Y = log( exp(Y0) + exp(Y1));
        estim_value{t}  = compute_value_function(Y, State_var);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %Policy Function (Investments)
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        %Store the optimal choice for each state (by P);
        Inv0 = I0(ind0);
        Inv1 = I1(ind1);

        Y0 = Inv0(:);
        estim_policy0{t} = compute_policy_function(Y0, State_var);
        Y1 = Inv1(:);
        estim_policy1{t} = compute_policy_function(Y1, State_var);


    end  %End "if last period or not";

end %End loop for "t";



sb.Hc_tp1_guess_prime = Hc_tp1_guess;
sb.estim_policy0 = estim_policy0;
sb.estim_policy1 = estim_policy1;
sb.estim_value0  = estim_value0;
sb.estim_value1  = estim_value1;
sb.Value0_actual = Value0_actual;
sb.Value1_actual = Value1_actual;
sb.distrib_mc_H_peers0 = distrib_mc_H_peers0;
sb.distrib_mc_H_peers1 = distrib_mc_H_peers1;